Unit wSimSpec;

Interface

uses Objects, WinTypes, WinProcs, WinDos, Wincrt, Strings,
     OStdDlgs, OWindows, ODialogs, OMemory,
     Cartocmds,
     tit_win,get_win,gmuf_win,aux_win,outl_win,coul_win;

Type
  pSpectre = ^Spectre;
  Spectre  = Array[1..751] of Double;
  pAcoWav  = ^AcoWav;
  Acowav   = Array[1..8,1..128] of Double;
  pRO_ARR = ^RO_ARR;
  RO_ARR = array[1..3,1..10] of double;
  pRM_ARR = ^RM_ARR;
  RM_ARR = array[1..34] of double;

  pEX_ARR = ^EX_ARR;
  EX_ARR = array[1..4,1..10] of double;
  pC_ARR = ^C_ARR;
  C_ARR = array[1..83] of double;
  pPH_ARR = ^PH_ARR;
  PH_ARR = array[1..10,1..83] of double;
  pEXT_ARR = ^EXT_ARR;
  EXT_ARR = Array[1..10] of Double;
  pCO3_ARR = ^CO3_ARR;
  CO3_ARR  = array[1..102] of Double;
  pCCH2O_ARR = ^CCH2O_ARR;
  CCH2O_ARR = array[1..15] of Double;

(* GLOBALS *)

var
 { Scattered direction }

  PHI,XMUS,XMUV,XMUP,XMUD,ADIF : Double;

 { Atmosphere model }

  UW,UO3       : Double;
  IDATM        : Integer;

const
  PI    : Double = 3.1415926;

(* Common FFU *)
var
  _S                  : pSpectre;
  _WLINF, _WLSUP      : Double;
(* Common DEL *)
const
      _SIGMA  : Double = 0.056032;
      _DELTA  : Double = 0.0139;
(* Common ATM *)
var
  Z,P,T,WH,WO         : pRM_ARR;
(* Common AER *)
var
  _EXT,OME,GASYM,PHASE : pEXT_ARR;
(* Common DISC *)
var
  ROATM,DTDIR,DTDIF,UTDIR,UTDIF,SPHAL : pRO_ARR;
  WLDIS,TRAYL                         : pEXT_ARR;
(*Common Aeroso *)
var
  PH_D,PH_W,PH_O,PH_S       : pPH_ARR;
  EX,SC,ASY                 : pEX_ARR;
  C_                        : pC_ARR;

var
(* Common ABS *)
  ACO          : pAcoWav;
  ACR_DICA1    : pAcoWav;
  ACR_DICA2    : pAcoWav;
  ACR_DICA3    : pAcoWav;
  ACR_OXYG3    : pAcoWav;
  ACR_OXYG4    : pAcoWav;
  ACR_OXYG5    : pAcoWav;
  ACR_OXYG6    : pAcoWav;
  ACR_OZON1    : pAcoWav;
  ACR_WAVA1    : pAcoWav;
  ACR_WAVA2    : pAcoWav;
  ACR_WAVA3    : pAcoWav;
  ACR_WAVA4    : pAcoWav;
  ACR_WAVA5    : pAcoWav;
  ACR_WAVA6    : pAcoWav;

var
  SI                  : pSpectre;
(* Spectral curves *)
  WLDISC              : pEXT_ARR;

var
   CO3    : pCO3_ARR;
   CCH2O  : pCCH2O_ARR;

FUNCTION   EXP_D(X:Double):Double;
procedure  ENVIRO(DIFR,DIFA,R:Double; var FRA,FAE,FR:Double);
procedure  INTERP (IAER:integer; WL,TAER55,XMUS,XMUV,XMUD:Double;
               var RORAYL,ROAERO,PHAA,PHAR,TSCA,
                   TRAY,TAER,DTOTT,UTOTT,
                   ASTOT,ASRAY,ASAER,
                   UTOTR,UTOTA,DTOTR,DTOTA:Double);
procedure  READ_WAVAP(var f_spec:text);
procedure  ABSTRA(IDATM:integer; WL,XMUS,XMUV: Double; var UW,UO3:Double;
              var DTWAVA,DTOZON,DTDICA,DTOXYG,
                  UTWAVA,UTOZON,UTDICA,UTOXYG,
                  TTWAVA,TTOZON,TTDICA,TTOXYG:Double);
procedure  DISCOM (IAER:integer;var XMUD,XMUS,XMUV,TAER55:Double);
procedure  READ_SOLIRR(var f_spec:text);
procedure  SOLIRR(WL:Double;var SWL:Double);
procedure  JINV;
procedure  ODA550(IAER:integer; V:Double; var TAER55:Double);
Procedure  READ_AEROSO(var f_spec:text);
Procedure  AEROSO (IAER:integer;var CO_IN; var XMUD:double);
Procedure  ATMOSPHERE_MODEL(iatmo:integer) ;
Function   ASIN(X:Double):Double;
Function   ACOS(X:Double):Double;
Function   POSSOL (MONTH,JDAY:integer;TU,XLON,XLAT:double;
                        var ASOL,PHI0:double):integer;
Function  POSLAN (MONTH,JDAY:integer;TU,XLON,XLAT:double;
                        var ASOL,PHI0,AVIS,PHIV:double):integer;
procedure READ_WLDISC(var f_spec:text);
procedure READ_CO3(var f_spec:text);
procedure READ_CCH2O(var f_spec:text);

Implementation

FUNCTION EXP_D(X:Double):Double;
(*
C********************************************************************
C
C	Special exponential function exp(x) which avoids arithmetic
C	overflow at large x or arithmetic underflow at small x
C
C********************************************************************
*)
var big,EX_ : Double;
BEGIN
	EX_ := 0.0;
	big := 85.0;
	if(X>big) then EX_ := EXP(big);
	IF(X>-big) and (X <= big) then EX_ := EXP(X);
	EXP_D := EX_;
END;

procedure ENVIRO(DIFR,DIFA,R:Double; var FRA,FAE,FR:Double);
(*C                                                                       *)
(*C     CALCULATION OF THE ENVIRONMENTAL FUNCTION FOR                     *)
(*C     RAYLEIGH AND AEROSOLS CONTRIBUTION.                               *)
(*C                                                                       *)
(*C     CAUTION !!!!                                                      *)
(*C                                                                       *)
(*C     THIS CALCULATION HAVE BEEN DONE FOR NADIR OBSERVATION             *)
(*C     BUT ONE MAY EXTRAPOLATE THE RESULT FOR ANGLE TILL 20 DEG.         *)
(*C                                                                       *)
BEGIN
      FAE:=1-0.375*EXP_D(-R*0.200)-0.625*EXP_D(-R*1.800);
      FRA:=1-0.930*EXP_D(-R*0.080)-0.070*EXP_D(-R*1.100);
      FR:=(FAE*DIFA+FRA*DIFR)/(DIFA+DIFR);
END;

procedure INTERP (IAER:integer; WL,TAER55,XMUS,XMUV,XMUD:Double;
               var RORAYL,ROAERO,PHAA,PHAR,TSCA,
                   TRAY,TAER,DTOTT,UTOTT,
                   ASTOT,ASRAY,ASAER,
                   UTOTR,UTOTA,DTOTR,DTOTA:Double);

(*C     THAT FOR THE ATMOSPHERE :                                         *)
(*C     THE REFLECTANCES                                                  *)
(*C                     RAYLEIGH                             := RORAYL     *)
(*C                     AEROSOLS                             := ROAERO     *)
(*C     THE DOWNWARD TRANSMITTANCES                                       *)
(*C                     RAYLEIGH                             := DTOTR      *)
(*C                     AEROSOLS                             := DTOTA      *)
(*C                     TOTAL                                := DTOTT      *)
(*C     THE UPWARD TRANSMITTANCES                                         *)
(*C                     RAYLEIGH                             := UTOTR      *)
(*C                     AEROSOLS                             := UTOTA      *)
(*C                     TOTAL                                := UTOTT      *)
(*C     THE SPHERICAL ALBEDOS                                             *)
(*C                     RAYLEIGH                             := ASRAY      *)
(*C                     AEROSOLS                             := ASAER      *)
(*C                     TOTAL                                := ASTOT      *)

var
  LINF,LL,LSUP : Integer;
  D2,ALPHAA,BETAA,ALPHAR,BETAR,ALPHAC,BETAC        : Double;
  DRINF,DRSUP,DTINF,DTSUP,DTOTC,DAINF,DASUP        : Double;
  URINF,URSUP,UTINF,UTSUP,UTOTC,UAINF,UASUP        : Double;
  ARINF,ARSUP,ATINF,ATSUP                          : Double;
  ASAEC,AAINF,AASUP,COEF_,WLINF_                   : Double;

label 1234,1235,1236,1237,1238,1239,1240;

BEGIN
      LINF:=1;
      for LL:=1 to 9 do begin
        IF(WL>=WLDIS^[LL]) AND (WL<=WLDIS^[LL+1]) then LINF:=LL;
      end;
      IF(WL>WLDIS^[1])then LINF:=9;
      LSUP:=LINF+1;
 
(*C     INTERPOLATION IN FUNCTION OF WAVELENGTH FOR SCATTERING            *)
(*C     ATMOSPHERIC FUNCTIONS FROM DISCRETE VALUES AT WLDIS               *)

      ALPHAA:=0.;
      BETAA:=0.;
      ALPHAR:=0.;
      BETAR:=0.;
      ALPHAC:=0.;
      BETAC:=0.;
      COEF_:=LN(WLDIS^[LSUP]/WLDIS^[LINF]);
      WLINF_:=WLDIS^[LINF];
(*C                                                                       *)
      IF(IAER=0) then GOTO 1240;
      ALPHAA:=LN(PHASE^[LSUP]/PHASE^[LINF])/COEF_;
      BETAA:=PHASE^[LINF]/(EXP(LN(WLINF_)*ALPHAA));
 1240:PHAA:=BETAA*(EXP(LN(WL)*ALPHAA));
      D2:=2.0+_DELTA;
      PHAR:=(2.0*(1.0-_DELTA)/D2)*0.75*(1.0+XMUD*XMUD)+3.0*_DELTA/D2;
      ALPHAR:=LN(ROATM^[1,LSUP]/ROATM^[1,LINF])/ COEF_;
      BETAR:=ROATM^[1,LINF]/(EXP(LN(WLINF_)*ALPHAR));
      IF(IAER=0) then GOTO 1234;
      ALPHAA:=LN(ROATM^[3,LSUP]/ROATM^[3,LINF])/COEF_;
      BETAA:=ROATM^[3,LINF]/(EXP(LN(WLINF_)*ALPHAA));
 1234:RORAYL:=BETAR*(EXP(LN(WL)*ALPHAR));
      ROAERO:=BETAA*(EXP(LN(WL)*ALPHAA));
      ALPHAR:=LN(TRAYL^[LSUP]/TRAYL^[LINF])/COEF_;
      BETAR:=TRAYL^[LINF]/(EXP(LN(WLINF_)*ALPHAR));
      TRAY:=BETAR*(EXP(LN(WL)*ALPHAR));
      IF(IAER=0) then GOTO 1235;
      ALPHAA:=LN(_EXT^[LSUP]*OME^[LSUP]/(_EXT^[LINF]*OME^[LINF]))/COEF_;
      BETAA:=_EXT^[LINF]*OME^[LINF]/(EXP(LN(WLINF_)*ALPHAA));
      TSCA:=TAER55*BETAA*(EXP(LN(WL)*ALPHAA))/_EXT^[4];
      ALPHAA:=LN(_EXT^[LSUP]/_EXT^[LINF])/COEF_;
      BETAA:=_EXT^[LINF]/(EXP(LN(WLINF_)*ALPHAA));
 1235:TAER:=TAER55*BETAA*(EXP(LN(WL)*ALPHAA))/_EXT^[4];
      DRINF:=DTDIF^[1,LINF]+DTDIR^[1,LINF];
      DRSUP:=DTDIF^[1,LSUP]+DTDIR^[1,LSUP];
      ALPHAR:=LN(DRSUP/DRINF)/COEF_;
      BETAR:=DRINF/(EXP(LN(WLINF_)*ALPHAR));
      DTOTR:=BETAR*(EXP(LN(WL)*ALPHAR));
      DTINF:=DTDIF^[2,LINF]+DTDIR^[2,LINF];
      DTSUP:=DTDIF^[2,LSUP]+DTDIR^[2,LSUP];
      ALPHAC:=LN((DTSUP*DRINF)/(DTINF*DRSUP))/COEF_;
      BETAC:=(DTINF/DRINF)/(EXP(LN(WLINF_)*ALPHAC));
      DTOTC:=BETAC*(EXP(LN(WL)*ALPHAC));
      DAINF:=DTDIF^[3,LINF]+DTDIR^[3,LINF];
      DASUP:=DTDIF^[3,LSUP]+DTDIR^[3,LSUP];
      IF(IAER=0) then GOTO 1236;
      ALPHAA:=LN(DASUP/DAINF)/COEF_;
      BETAA:=DAINF/(EXP(LN(WLINF_)*ALPHAA));
 1236:DTOTA:=BETAA*(EXP(LN(WL)*ALPHAA));
      DTOTT:=DTOTC*DTOTR;
      URINF:=UTDIF^[1,LINF]+UTDIR^[1,LINF];
      URSUP:=UTDIF^[1,LSUP]+UTDIR^[1,LSUP];
      ALPHAR:=LN(URSUP/URINF)/ COEF_;
      BETAR:=URINF/(EXP(LN(WLINF_)*ALPHAR));
      UTOTR:=BETAR*(EXP(LN(WL)*ALPHAR));
      UTINF:=UTDIF^[2,LINF]+UTDIR^[2,LINF];
      UTSUP:=UTDIF^[2,LSUP]+UTDIR^[2,LSUP];
      ALPHAC:=LN((UTSUP*URINF)/(UTINF*URSUP))/ COEF_;
      BETAC:=(UTINF/URINF)/(EXP(LN(WLINF_)*ALPHAC));
      UTOTC:=BETAC*(EXP(LN(WL)*ALPHAC));
      UAINF:=UTDIF^[3,LINF]+UTDIR^[3,LINF];
      UASUP:=UTDIF^[3,LSUP]+UTDIR^[3,LSUP];
      IF(IAER=0) then GOTO 1237;
      ALPHAA:=LN(UASUP/UAINF)/ COEF_;
      BETAA:=UAINF/(EXP(LN(WLINF_)*ALPHAA));
 1237:UTOTA:=BETAA*(EXP(LN(WL)*ALPHAA));
      UTOTT:=UTOTC*UTOTR;
      ARINF:=SPHAL^[1,LINF];
      ARSUP:=SPHAL^[1,LSUP];
      ALPHAR:=LN(ARSUP/ARINF)/ COEF_;
      BETAR:=ARINF/(EXP(LN(WLINF_)*ALPHAR));
      ASRAY:=BETAR*(EXP(LN(WL)*ALPHAR));
      ATINF:=SPHAL^[2,LINF];
      ATSUP:=SPHAL^[2,LSUP];
      ALPHAC:=0.0;
      BETAC:=0.0;
      IF(IAER=0) THEN GOTO 1238;
      ALPHAC:=LN((ATSUP-ARSUP)/(ATINF-ARINF))/COEF_;
      BETAC:=(ATINF-ARINF)/(EXP(LN(WLINF_)*ALPHAC));
      ASAEC:=BETAC*(EXP(LN(WL)*ALPHAC));
 1238:
      AAINF:=SPHAL^[3,LINF];
      AASUP:=SPHAL^[3,LSUP];
      IF(IAER=0) then GOTO 1239;
      ALPHAA:=LN(AASUP/AAINF)/COEF_;
      BETAA:=AAINF/(EXP(LN(WLINF_)*ALPHAA));
 1239:ASAER:=BETAA*(EXP(LN(WL)*ALPHAA));
      ASTOT:=ASAEC+ASRAY;

 END;

procedure READ_WAVAP_ONE(var f_spec:text;ACR:pAcoWav);
var nbrl,i,j,k : integer;
  erreur       : Integer;
begin

  j := 0;
  repeat
     readln(f_spec,S);
     if copy(S,length(S)-2,3) = '*0.' then begin
       S := copy(S,1,length(S)-3);
       val(S, nbrl, erreur); nbrl := nbrl div 8;
       for i := 1 to nbrl do begin
         j:=j+1;
 (*        writeln(j);    *)
         for k := 1 to 8 do ACR^[k,j]:= 0;
       end;
     end else begin
       j:=j+1;
       for k := 1 to 5 do begin
          _nextpar(strw,S); val(strw, ACR^[k,j], erreur);
       end;
       readln(f_spec,S); for k := 1 to 3 do begin
          _nextpar(strw,S); val(strw, ACR^[k+5,j], erreur);
       end;
(*       writeln(j,' ',ACR^[6,j]);  *)
     end;
  until j = 128;

(*  writeln(ACR^[6,j],' ',ACR^[7,j],' ',ACR^[8,j]);  *)

end;

procedure READ_WAVAP(var f_spec:text);
begin

 READ_WAVAP_ONE(f_spec, ACR_DICA1);
 READ_WAVAP_ONE(f_spec, ACR_DICA2);
 READ_WAVAP_ONE(f_spec, ACR_DICA3);
 READ_WAVAP_ONE(f_spec, ACR_OXYG3);
 READ_WAVAP_ONE(f_spec, ACR_OXYG4);
 READ_WAVAP_ONE(f_spec, ACR_OXYG5);
 READ_WAVAP_ONE(f_spec, ACR_OXYG6);
 READ_WAVAP_ONE(f_spec, ACR_OZON1);
 READ_WAVAP_ONE(f_spec, ACR_WAVA1);
 READ_WAVAP_ONE(f_spec, ACR_WAVA2);
 READ_WAVAP_ONE(f_spec, ACR_WAVA3);
 READ_WAVAP_ONE(f_spec, ACR_WAVA4);
 READ_WAVAP_ONE(f_spec, ACR_WAVA5);
 READ_WAVAP_ONE(f_spec, ACR_WAVA6);

end;

procedure ABSTRA(IDATM:integer; WL,XMUS,XMUV: Double; var UW,UO3:Double;
              var DTWAVA,DTOZON,DTDICA,DTOXYG,
                  UTWAVA,UTOZON,UTDICA,UTOXYG,
                  TTWAVA,TTOZON,TTDICA,TTOXYG:Double);

(*C     TRANSMITTANCE CALCULATION FOR OZONE  WATER VAPOR                  *)
(*C     CARBON DIOXYDE AND OXYGEN.                                        *)
 
(*C     DOWNWARD ABSORPTION WATER VAPOR DTWAVA                            *)
(*C     DOWNWARD ABSORPTION OZONE       DTOZON                            *)
(*C     DOWNWARD ABSORPTION CARBON DIOX DTDICA                            *)
(*C     DOWNWARD ABSORPTION OXYGEN      DTOXYG                            *)
(*C     UPWARD ABSORPTION WATER VAPOR   UTWAVA                            *)
(*C     UPWARD ABSORPTION OZONE         UTOZON                            *)
(*C     UPWARD ABSORPTION CARBON DIOX   UTDICA                            *)
(*C     UPWARD ABSORPTION OXYGEN        UTOXYG                            *)
(*C     TOTAL(ON THE TWO PATHS ) ABSORPTION WATER VAPOR TTWAVA            *)
(*C     TOTAL(ON THE TWO PATHS ) ABSORPTION OZONE       TTOZON            *)
(*C     TOTAL(ON THE TWO PATHS ) ABSORPTION CARBON DIOX TTDICA            *)
(*C     TOTAL(ON THE TWO PATHS ) ABSORPTION OXYGEN      TTOXYG            *)

procedure DICA1;
(*C                                                                       *)
(*C    CARBON DIOXYDE : ( 2500 - 5040 CM-1 )                              *)
(*C    ACR_DICA1                                                          *)
var k,j : integer;
begin
  for k := 1 to 8 do
    for j := 1 to 128 do ACO^[k,j] := ACR_DICA1^[k,j];
end;

procedure DICA2;
(*C                                                                       *)
(*C    CARBON DIOXYDE : ( 5060 - 7600 CM-1 )                              *)
(*C                                                                       *)
var k,j : integer;
begin
  for k := 1 to 8 do
    for j := 1 to 128 do ACO^[k,j] := ACR_DICA2^[k,j];
end;

procedure DICA3;
(*C                                                                       *)
(*C    CARBON DIOXYDE : ( 7620 - 9620 CM-1 )                              *)
(*C                                                                       *)
var k,j : integer;
begin
  for k := 1 to 8 do
    for j := 1 to 128 do ACO^[k,j] := ACR_DICA3^[k,j];
end;

procedure OXYG3;
(*C                                                                       *)
(*C     OXYGEN : ( 7620 - 10160 CM-1 )                                    *)
(*C                                                                       *)
var k,j : integer;
begin
  for k := 1 to 8 do
    for j := 1 to 128 do ACO^[k,j] := ACR_OXYG3^[k,j];
end;

procedure OXYG4;
(*C                                                                       *)
(*C     OXYGEN : ( 10180 - 12720 CM-1 )                                   *)
(*C                                                                       *)
var k,j : integer;
begin
  for k := 1 to 8 do
    for j := 1 to 128 do ACO^[k,j] := ACR_OXYG4^[k,j];
end;

procedure OXYG5;
(*C                                                                       *)
(*C     OXYGEN : ( 12740 - 15280 CM-1 )                                   *)
(*C                                                                       *)
var k,j : integer;
begin
  for k := 1 to 8 do
    for j := 1 to 128 do ACO^[k,j] := ACR_OXYG5^[k,j];
end;

procedure OXYG6;
(*C                                                                       *)
(*C     OXYGEN : ( 15300 - 15900 CM-1 )                                   *)
(*C                                                                       *)
var k,j : integer;
begin
  for k := 1 to 8 do
    for j := 1 to 128 do ACO^[k,j] := ACR_OXYG6^[k,j];
end;

procedure OZON1;
(*C                                                                       *)
(*C     OZONE :  ( 2500 - 3040 CM-1 )                                     *)
(*C                                                                       *)
var k,j : integer;
begin
  for k := 1 to 8 do
    for j := 1 to 128 do ACO^[k,j] := ACR_OZON1^[k,j];
end;

procedure WAVA1;
(*C                                                                       *)
(*C     WATER VAPOR (2500 - 5040 CM-1)                                    *)
(*C                                                                       *)
var k,j : integer;
begin
  for k := 1 to 8 do
    for j := 1 to 128 do ACO^[k,j] := ACR_WAVA1^[k,j];
end;

procedure WAVA2;
(*C                                                                       *)
(*C     WATER VAPOR (5060 - A600 CM-1)                                    *)
(*C                                                                       *)
var k,j : integer;
begin
  for k := 1 to 8 do
    for j := 1 to 128 do ACO^[k,j] := ACR_WAVA2^[k,j];
end;

procedure WAVA3;
(*C                                                                       *)
(*C     WATER VAPOR (7620 - 10160 CM-1)                                   *)
(*C                                                                       *)
var k,j : integer;
begin
  for k := 1 to 8 do
    for j := 1 to 128 do ACO^[k,j] := ACR_WAVA3^[k,j];
end;

procedure WAVA4;
(*C                                                                       *)
(*C     WATER VAPOR (10180 - 12720 CM-1)                                  *)
(*C                                                                       *)
var k,j : integer;
begin
  for k := 1 to 8 do
    for j := 1 to 128 do ACO^[k,j] := ACR_WAVA4^[k,j];
end;

procedure WAVA5;
(*C                                                                       *)
(*C     WATER VAPOR (12740 - 15280 CM-1)                                  *)
(*C                                                                       *)
var k,j : integer;
begin
  for k := 1 to 8 do
    for j := 1 to 128 do ACO^[k,j] := ACR_WAVA5^[k,j];
end;

procedure WAVA6;
(*C                                                                       *)
(*C     WATER VAPOR (15300 - 17860 CM-1)                                  *)
(*C                                                                       *)
var k,j : integer;
begin
  for k := 1 to 8 do
    for j := 1 to 128 do ACO^[k,j] := ACR_WAVA6^[k,j];
end;

const
   IVLI : array[1..6] of integer = (2500, 5060, 7620, 10180, 12740, 15300);
{
   CO3  : array[1..102] of Double =  (
      4.50E-03,  8.00E-03,  1.07E-02,  1.10E-02,  1.27E-02,  1.71E-02,
      2.00E-02,  2.45E-02,  3.07E-02,  3.84E-02,  4.78E-02,  5.67E-02,
      6.54E-02,  7.62E-02,  9.15E-02,  1.00E-01,  1.09E-01,  1.20E-01,
      1.28E-01,  1.12E-01,  1.11E-01,  1.16E-01,  1.19E-01,  1.13E-01,
      1.03E-01,  9.24E-02,  8.28E-02,  7.57E-02,  7.07E-02,  6.58E-02,
      5.56E-02,  4.77E-02,  4.06E-02,  3.87E-02,  3.82E-02,  2.94E-02,
      2.09E-02,  1.80E-02,  1.91E-02,  1.66E-02,  1.17E-02,  7.70E-03,
      6.10E-03,  8.50E-03,  6.10E-03,  3.70E-03,  3.20E-03,  3.10E-03,
      2.55E-03,  1.98E-03,  1.40E-03,  8.25E-04,  2.50E-04,  0.,
      0.,        0.,        5.65E-04,  2.04E-03,  7.35E-03,  2.03E-02,
      4.98E-02,  1.18E-01,  2.46E-01,  5.18E-01,  1.02E+00,  1.95E+00,
      3.79E+00,  6.65E+00,  1.24E+01,  2.20E+01,  3.67E+01,  5.95E+01,
      8.50E+01,  1.26E+02,  1.68E+02,  2.06E+02,  2.42E+02,  2.71E+02,
      2.91E+02,  3.02E+02,  3.03E+02,  2.94E+02,  2.77E+02,  2.54E+02,
      2.26E+02,  1.96E+02,  1.68E+02,  1.44E+02,  1.17E+02,  9.75E+01,
      7.65E+01,  6.04E+01,  4.62E+01,  3.46E+01,  2.52E+01,  2.00E+01,
      1.57E+01,  1.20E+01,  1.00E+01,  8.80E+00,  8.30E+00,  8.60E+00);
const
   CCH2O : array[1..15] of Double =  (
     10.00, 0.19, 0.15, 0.12, 0.10, 0.09, 0.10, 0.12, 0.15, 0.17, 0.20,
     0.24, 20.28, 0.33, 0.00);
}

type
  pTNU_ARR = ^TNU_ARR;
  TNU_ARR = array[1..4,1..3] of double;
var
  TNU            : pTNU_ARR;
  A              : array[1..8] of double;
  RM,R2,R3,TP    : pRM_ARR;
  RAT            : array[1..4] of double;
  DTCONT,UTCONT,TTCONT,PSI,DS2 : Double;
  V,P0,G,T0,DS,TE,ROAIR,AIR,ROCO2,RMO2,RMO3,UWUS,UO3US : Double;
  IV,ID,IDGAZ,N,NH,INU,II : integer;
  UU,U,UP,UUD,UUT,UUU,UD,UT,UPD,UPT,UDP,UPDP,UDTP,UPDTP : Double;
  XD,XI,XH,AH2O,AKO3,TEST1,TEST2,TEST3,UDT,UPDT,TN,ATEST : Double;
  TT,Y,UTT,UPTT,PHI : Double;
  i,j,k : integer;

label 40,45,130,135,140,145,146,147,148,270,271,999;

begin              (* procedure ABSTRA *)

  getmem(TNU,sizeOf(TNU_ARR));
  getmem(RM,sizeOf(RM_ARR));
  getmem(R2,sizeOf(RM_ARR));
  getmem(R3,sizeOf(RM_ARR));
  getmem(TP,sizeOf(RM_ARR));

      DTWAVA:=1.;
      UTWAVA:=1.;
      TTWAVA:=1.;
      DTCONT:=1.;
      UTCONT:=1.;
      TTCONT:=1.;
      DTOZON:=1.;
      UTOZON:=1.;
      TTOZON:=1.;
      DTDICA:=1.;
      UTDICA:=1.;
      TTDICA:=1.;
      DTOXYG:=1.;
      UTOXYG:=1.;
      TTOXYG:=1.;
      for I:=1 to 4 do RAT[I]:=1.0;
      IF (IDATM = 0) then goto 999;
      IF (XMUS=0.0) OR (XMUV=0.0) then begin
        Loadstring(HInstance, stW_SpecZ90+Lang, MessageText, 79);
        Loadstring(HInstance, stW_SpecTIT+Lang, TitleText, 79);
        MessageBox(0, MessageText, TitleText, mb_IconStop);
        GOTO 999;
      end;

(*C     CONSTANTS DETERMINATION                                           *)
 
      P0:=1013.25;
      G:=98.1;
      T0:=250.;
 
(*C     VOLUMIC MASS IN KILOGRAMMES PER M3                                *)
 
      DS:=0.0;
      TE:=0.0;
      ROAIR:=0.0;
      AIR:=0.028964/0.0224;
      ROCO2:=0.044/0.0224;
      RMO2:=0.032/0.0224;
      RMO3:=0.048/0.0224;
      UWUS:=1.424;
      UO3US:=0.344;

(*      writeln('ABSTRA - '+'IDATM : ',IDATM,' WL ',WL:7:4);  *)
      IF(IDATM=8) then begin
        RAT[1]:=UW/UWUS;
        RAT[2]:=1.0;
        RAT[3]:=1.0;
        RAT[4]:=UO3/UO3US;
      end;
      V:=1.0E+04/WL;
      IV:=round(V/5.0);
      IV:=IV*5;
      ID:=((IV-2500) div 20) div 128+1;
(*      writeln('ABSTRA - '+'IV : ',IV,'ID : ',ID);   *)

 For IDGAZ := 1 to 4 do begin

(*      writeln('ABSTRA - IDGAZ ',IDGAZ);    *)
      if not(ID in [1..6]) then goto 270;
      case ID of
        1: begin   (* 101 *)
             case IDGAZ of
               1: WAVA1;
               2: DICA1;
               3: goto 270;
               4: OZON1;
             end;
           end;
        2: begin   (* 102 *)
             case IDGAZ of
               1: WAVA2;
               2: DICA2;
               3: goto 270;
               4: goto 270;
             end;
           end;
        3: begin   (* 103 *)
             case IDGAZ of
               1: WAVA3;
               2: DICA3;
               3: OXYG3;
               4: goto 270;
             end;
           end;
        4: begin   (* 104 *)
             case IDGAZ of
               1: WAVA4;
               2: goto 270;
               3: OXYG4;
               4: goto 270;
             end;
           end;
        5: begin   (* 105 *)
             case IDGAZ of
               1: WAVA5;
               2: goto 270;
               3: OXYG5;
               4: goto 270;
             end;
           end;
        6: begin   (* 106 *)
             case IDGAZ of
               1: WAVA6;
               2: goto 270;
               3: OXYG6;
               4: goto 270;
             end;
           end;
      end;
      INU:=(IV-IVLI[ID]) div 20+1;
      for II:=1 to 8 do A[II]:=ACO^[II,INU];
      GOTO 271;
270:  for II := 1 to 8 do A[II]:= 0.0;

(*C     MIXING RATIO CALCULATION FOR EACH GASEOUS CONSTITUENTS            *)

271:
(*
      writeln('ABSTRA - A : ',A[1]:7:4,' ',A[2]:7:4,' ',A[3]:7:4,' ',A[4]:7:4);
      writeln('ABSTRA - A : ',A[5]:7:4,' ',A[6]:7:4,' ',A[7]:7:4,' ',A[8]:7:4);
*)
     for K:=1 to 33 do begin
        ROAIR:=AIR*273.16*P^[K]/(1013.25*T^[K]);
        IF(IDGAZ=1)then RM^[K]:=WH^[K]/(ROAIR*1000.0);
        IF(IDGAZ=2)then RM^[K]:=3.3E-04*ROCO2/AIR;
        IF(IDGAZ=3)then RM^[K]:=0.20947*RMO2/AIR;
        IF(IDGAZ=4)then RM^[K]:=WO^[K]/(ROAIR*1000.0);
        TP^[K]:=(T^[K]+T^[K+1])/2.;
        TE:=TP^[K]-T0;
        PHI:=EXP_D(A[3]*TE+A[4]*(TE*TE));
        PSI:=EXP_D(A[5]*TE+A[6]*(TE*TE));
        R2^[K]:=RM^[K]*PHI;
        R3^[K]:=RM^[K]*PSI;
      end;
(*      writeln('ABSTRA - Initialisation done'); *)

      UU:=0.;
      U:=0.;
      UP:=0.;
      UUD:=0.;
      UUT:=0.;
      UUU:=0.;
      UD:=0.;
      UT:=0.;
      UPD:=0.;
      UPT:=0.;
      UDP:=0.;
      UPDP:=0.;
      UDTP:=0.;
      UPDTP:=0.;

(*C     PRESSURE SCALE INTEGRATION SIGMA:=P/P0                             *)
 
      for K:= 2 to 33 do begin
        DS:=(P^[K-1]-P^[K])/P^[1];
        DS2:=(P^[K-1]*P^[K-1]-P^[K]*P^[K])/(2.0*P^[1]*P0);
        UU:=UU+((RM^[K]+RM^[K-1])/2.0)*DS*RAT[IDGAZ];
        U:=U+((R2^[K]+R2^[K-1])/2.0)*DS*RAT[IDGAZ];
        UP:=UP+((R3^[K]+R3^[K-1])/2.0)*DS2*RAT[IDGAZ];
      end;
(*      writeln('ABSTRA - Pressure scale integration done');  *)

      UU:=UU*P^[1]*100.0/G;
      U:=U*P^[1]*100.0/G;
      UP:=UP*P^[1]*100.0/G;
 
      IF(IDGAZ=4) then UU:=1000*UU/RMO3;
      IF(IDGAZ=2) then UU:=1000*UU/ROCO2;
 
(*C     DOWNWARD PATH                                                     *)
 
      UUD:=UU/XMUS;
 
(*C     UPWARD PATH                                                       *)
 
      UUU:=UU/XMUV;
 
(*C     TOTAL(DOWN+UP) PATH                                               *)

      UUT:=UU*(1./XMUS+1./XMUV);
      IF(IDGAZ=1) then GOTO 146;
      IF(IDGAZ=2) AND (IV>9620)then GOTO 147;
      IF(IDGAZ=3) AND (IV>15920)then GOTO 147;
      IF(IDGAZ=4) then GOTO 146;
      GOTO 145;
  146:XI:=0.;
      N:=0;
      XD:=0.;
      IF(IV<2350) OR (IV>3000) then GOTO 148;
      XI:=(V-2350.0)/50.0+1.0;
      NH:=round(XI+1.001);
      XH:=XI-NH;
      AH2O:=CCH2O^[NH]+XH*(CCH2O^[NH]-CCH2O^[NH-1]);
      DTCONT:=EXP_D(-AH2O*UUD);
      UTCONT:=EXP_D(-AH2O*UUU);
      TTCONT:=EXP_D(-AH2O*UUT);
  148:IF (IDGAZ=1) then GOTO 145;
      IF(IV<13000) then GOTO 145;
      IF(IV<=23400) then GOTO 130;
      IF(IV>=27500) then GOTO 135;
      TNU^[4,1]:=1.0;
      TNU^[4,2]:=1.0;
      TNU^[4,3]:=1.0;
      GOTO 45;
  130:XI:=(V-13000.0)/200.0+1.0;
      GOTO 140;
  135:XI:=(V-27500.0)/500.0+57.0;
  140:N:=round(XI+1.001);
      XD:=XI-N;
      AKO3:=CO3^[N]+XD*(CO3^[N]-CO3^[N-1]);
      TEST1:=AKO3*UUD;
      TEST2:=AKO3*UUU;
      TEST3:=AKO3*UUT;

(*C     EXPONENTIAL OVERFLOW TEST                                         *)
 
      IF(TEST1>86.0) then TEST1:=86.0;
      IF(TEST2>86.0) then TEST2:=86.0;
      IF(TEST3>86.0) then TEST3:=86.0;
 
      TNU^[4,1]:=EXP_D(-TEST1);
      TNU^[4,2]:=EXP_D(-TEST2);
      TNU^[4,3]:=EXP_D(-TEST3);
      GOTO 40;

  145:IF(IDGAZ=1) AND (IV>17860) then GOTO 147;
      IF(IDGAZ=4) AND (IV>3020) then GOTO 147;
(*C                                                                       *)
(*C     DOWNWARD PATH                                                     *)
(*C                                                                       *)
      UD:=U/XMUS;
      UPD:=UP/XMUS;
      UDT:=UD;
      IF(UD=0) AND (UPD=0.0) then UDT:=1.;
      TN:=A[2]*UPD/(2*UDT);
      ATEST:=A[2];
      IF (A[2]=0.0) AND (A[1]=0.0)then ATEST:=1.;
      UPDT:=UPD;
      IF(UD=0.0) AND (UPD=0.0)then UPDT:=1.;
      TT:=1+4*(A[1]/ATEST)*((UD*UD)/UPDT);
      Y:=-TN*(SQRT(TT)-1);
      IF(IDGAZ=1) then Y:=-A[1]*UD/SQRT(1+(A[1]/ATEST)*(UD*UD/UPDT));
      TNU^[IDGAZ,1]:=EXP_D(Y);
(*C                                                                       *)
(*C     UPWARD PATH                                                       *)
(*C                                                                       *)
      UDP:=U/XMUV;
      UPDP:=UP/XMUV;
      UDTP:=UDP;
      IF(UDP=0.0) AND (UPDP=0.0) then UDTP:=1.0;
      TN:=A[2]*UPDP/(2*UDTP);
      ATEST:=A[2];
      IF (A[2]=0.0) AND (A[1]=0.0) then ATEST:=1.0;
      UPDTP:=UPDP;
      IF(UDP=0.0) AND (UPDP=0.0)then UPDTP:=1.0;
      TT:=1+4*(A[1]/ATEST)*((UDP*UDP)/UPDTP);
      Y:=-TN*(SQRT(TT)-1);
      IF(IDGAZ=1) then Y:=-A[1]*UDP/SQRT(1+(A[1]/ATEST)*(UDP*UDP/UPDTP));
      TNU^[IDGAZ,2]:=EXP_D(Y);
(*C                                                                       *)
(*C     TOTAL(DOWN+UP) PATH                                               *)
(*C                                                                       *)
      UT:=U*(1.0/XMUS+1.0/XMUV);
      UPT:=UP*(1.0/XMUS+1.0/XMUV);
      UTT:=UT;
      IF(UT=0.0) AND (UPT=0.0)then UTT:=1.;
      TN:=A[2]*UPT/(2*UTT);
      UPTT:=UPT;
      IF(UT=0.0) AND (UPT=0.0) then UPTT:=1.;
      TT:=1+4*(A[1]/ATEST)*((UT*UT)/UPTT);
      Y:=-TN*(SQRT(TT)-1);
      IF(IDGAZ=1) then Y:=-A[1]*UT/SQRT(1+(A[1]/ATEST)*(UT*UT/UPTT));
      TNU^[IDGAZ,3]:=EXP_D(Y);
      GOTO 40;
  147:TNU^[IDGAZ,1]:=1.;
      TNU^[IDGAZ,2]:=1.;
      TNU^[IDGAZ,3]:=1.;
40:
 end;
   45:DTWAVA:=TNU^[1,1]*DTCONT;
      UTWAVA:=TNU^[1,2]*UTCONT;
      TTWAVA:=TNU^[1,3]*TTCONT;
      DTDICA:=TNU^[2,1];
      UTDICA:=TNU^[2,2];
      TTDICA:=TNU^[2,3];
      DTOXYG:=TNU^[3,1];
      UTOXYG:=TNU^[3,2];
      TTOXYG:=TNU^[3,3];
      DTOZON:=TNU^[4,1];
      UTOZON:=TNU^[4,2];
      TTOZON:=TNU^[4,3];

999:end;

procedure DISCOM (IAER:integer;var XMUD,XMUS,XMUV,TAER55:Double);
(*C     COMPUTATION OF ALL SCATTERING PARAMETERS AT WAVELENGTH            *)
(*C     DISCRETE VALUES SO WE                                             *)
(*C     CAN INTERPOLATE AT ANY WAVELENGTH                                 *)

procedure ODRAYL(WL:Double;var TRAY:Double);
(*C     MOLECULAR OPTICAL DEPTH                                           *)
Type
  pAT_ARR = ^AT_ARR;
  AT_ARR = array[1..70] of double;
  pDP_ARR = ^DP_ARR;
  DP_ARR = array[1..2,1..34] of double;

var
      A1, A2, A3, A4, AWL, AN, ANS, FNZ, FBETA, ANZ, BETAZ, A, AK : Double;
      AT, AP : pAT_ARR;
      DP, DT : pDP_ARR;
      DPZ, DTZ : Double;
      i,j,k,l : integer;
const
      PI  : Double = 3.1415926;
BEGIN

      getmem(AT, sizeof(AT_ARR));
      getmem(AP, sizeof(AT_ARR));
      getmem(DP, sizeof(DP_ARR));
      getmem(DT, sizeof(DP_ARR));

      AK:=1.0/WL;
      AWL:=WL;
(*C     AIR REFRACTION INDEX                                              *)
      A1:=146.0-AK*AK;
      A2:=41.0-AK*AK;
      A3:=3013189.3/A1;
      A4:=26113.82/A2;
      AN:=(6593.51+A3+A4)*1.0E-08;
      AN:=AN+1.E+00;
      A:=(8.0*PI*PI*PI/3.0)*
         ((AN*AN-1.0)*(AN*AN-1.0))*(6.0+3.0*_DELTA)/(6.0-7.0*_DELTA);
      FNZ:=1.0032882;
      FBETA:=A*FNZ*1E+01/(AWL*AWL*AWL*AWL);
      TRAY:=0.0;
      DP^[1,34]:=0.0;
      DP^[2,34]:=0.0;
      for K:=2 to 33 do begin
        DP^[1,K]:=EXP_D((LN(P^[K])+LN(P^[K-1]))/2.0);
        DP^[2,K]:=P^[K];
      end;
      for K:=2 to 34 do begin
        DT^[1,K]:=(T^[K]+T^[K-1])/2.;
        DT^[2,K]:=T^[K];
      end;
      L:=1;
      AP^[1]:=P^[1];
      AT^[1]:=T^[1];
      for K:=2 to 33 do begin
        for J:=1 to 2 do begin
          L:=L+1;
          AP^[L]:=DP^[J,K];
          AT^[L]:=DT^[J,K];
        end;
      end;

      for K:=1 to 51 do begin
        DPZ:=EXP_D((LN(AP^[K+1])+LN(AP^[K]))/2.0);
        DTZ:=(AT^[K+1]+AT^[K])/2.;
        ANZ:=DPZ/DTZ;
        BETAZ:=FBETA*ANZ;
        TRAY:=TRAY+BETAZ*0.5;
      end;
      for K:=52 to 61 do begin
        DPZ:=EXP_D((LN(AP^[K+1])+LN(AP^[K]))/2.0);
        DTZ:=(AT^[K+1]+AT^[K])/2.0;
        ANZ:=DPZ/DTZ;
        BETAZ:=FBETA*ANZ;
        TRAY:=TRAY+BETAZ*2.5;
      end;
      for K:=62 to 64 do begin
        DPZ:=EXP_D((LN(AP^[K+1])+LN(AP^[K]))/2.0);
        DTZ:=(AT^[K+1]+AT^[K])/2.0;
        ANZ:=DPZ/DTZ;
        BETAZ:=FBETA*ANZ;
        TRAY:=TRAY+BETAZ*5.0;
      end;
      ANZ:=P^[1]/T^[1];
      BETAZ:=FBETA*ANZ;
      TRAY:=TRAY+BETAZ*0.5;
(*
      str(WL:7:4,strw); _write(' WL : '+Strw);
      str(TRAY:7:4,strw); _writeln(' TRAY : '+Strw);
*)
      freemem(AT, sizeof(AT_ARR));
      freemem(AP, sizeof(AT_ARR));
      freemem(DP, sizeof(DP_ARR));
      freemem(DT, sizeof(DP_ARR));

END;

procedure ATMREF(IAER,XMUD,XMUS,XMUV,TRAYL_,TAER,GA,PIZA,PFA:Double;
                  var RORAYL,ROAERO:Double);
var D2,PHAR : Double;
    AK2,AK,E,F,DP,D,B,DEL,WW,SS : Double;
    Q1,Q2,Q3,C1,C2,CP1,CP2,Z,X,Y,AA1,AA2,AA3 :Double;

label 999;

BEGIN
(*C                                                                       *)
(*C     ATMOSPHERIC REFLECTANCES                                          *)
(*C                                                                       *)
      RORAYL:=0.0;
      ROAERO:=0.0;
(*C                                                                       *)
(*C     RAYLEIGH PHASE FUNCTION                                           *)
(*C                                                                       *)
      D2:=2.0+_DELTA;
      PHAR:=(2.0*(1.0-_DELTA)/D2)*0.75*(1.0+XMUD*XMUD)+(3.0*_DELTA)/D2;
(*C                                                                       *)
(*C     RAYLEIGH REFLECTANCE                                              *)
(*C                                                                       *)
      RORAYL:=TRAYL_*PHAR/(4.0*XMUS*XMUV);
      IF (IAER=0) then goto 999;
(*C                                                                       *)
(*C     AEROSOL REFLECTANCE                                               *)
(*C                                                                       *)
      AK2:=(1.0-PIZA)*(3.0-PIZA*3*GA);
      AK:=SQRT(AK2);
      E:=-3.0*XMUS*XMUS*PIZA/(4.0*(1.-AK2*XMUS*XMUS));
      F:=-(1.0-PIZA)*3*GA*XMUS*XMUS*PIZA/(4.0*(1.0-AK2*XMUS*XMUS));
      DP:=E/(3.0*XMUS)+XMUS*F;
      D:=E+F;
      B:=2*AK/(3.0-PIZA*3*GA);
      DEL:=EXP_D(AK*TAER)*(1.0+B)*(1.0+B)-EXP_D(-AK*TAER)*(1-B)*(1-B);
      WW:=PIZA/4.0;
      SS:=XMUS/(1.-AK2*XMUS*XMUS);
      Q1:=2.0+3.0*XMUS+(1.-PIZA)*3*GA*XMUS*(1.0+2.0*XMUS);
      Q2:=2.0-3.0*XMUS-(1.-PIZA)*3*GA*XMUS*(1.0-2.0*XMUS);
      Q3:=Q2*EXP_D(-TAER/XMUS);
      C1:=((WW*SS)/DEL)*(Q1*EXP(AK*TAER)*(1+B)+Q3*(1.0-B));
      C2:=-((WW*SS)/DEL)*(Q1*EXP(-AK*TAER)*(1-B)+Q3*(1.0+B));
      CP1:=C1*AK/(3.0-PIZA*3*GA);
      CP2:=-C2*AK/(3.0-PIZA*3*GA);
      Z:=D-PIZA*3*GA*XMUV*DP+PIZA*PFA/4.0;
      X:=C1-PIZA*3*GA*XMUV*CP1;
      Y:=C2-PIZA*3*GA*XMUV*CP2;
      AA1:=XMUV/(1.0+AK*XMUV);
      AA2:=XMUV/(1.0-AK*XMUV);
      AA3:=XMUS*XMUV/(XMUS+XMUV);
      ROAERO:=(1.0/(XMUS*XMUV))*(X*AA1*(1.-EXP_D(-TAER/AA1))
       +Y*AA2*(1-EXP_D(-TAER/AA2))+Z*AA3*(1.-EXP_D(-TAER/AA3)));

999:END;

FUNCTION TRANS(XMU,B1,PIZ,TAUT:Double):Double;
(*C                                                                       *)
(*C     COMPUTATIONS OF THE TOTAL TRANSMITTANCE                           *)
(*C                                                                       *)
(*      COMMON /TRANSM/B1 PIZ TAUT;    *)
var
    Z1,Z2,Z3,Z4,O2,AK,G,R6,E,H,R7,A,B,S7 : Double;
    R8,R9,Q7,X1,X2,AL4 : Double;
BEGIN
      IF ((1.0-PIZ)<1.0E-07) then PIZ:=1.0-1.0E-07;
      Z1:=2.0*(1.0-PIZ*(1.0-((3.0-B1)/8.0)))/
              (1.0-PIZ*(B1/3.0)*(B1/3.0));
      Z2:=2.0*PIZ*((3.0-B1)/8.0)/
              (1.0-PIZ*(B1/3.0)*(B1/3.0));
      O2:=0.5-XMU*(B1-3.0*(B1/3.0)*(B1/3.0))/
              4.0/(1.0-(B1/3.0)*(B1/3.0));
      Z3:=(1.0-(B1/3.0)*(B1/3.0))*PIZ*O2/
              (1.0-PIZ*(B1/3.0)*(B1/3.0));
      Z4:=(1.0-(B1/3.0)*(B1/3.0))*PIZ*(1.0-O2)/
              (1.0-PIZ*(B1/3.0)*(B1/3.0));
      AK:=SQRT(Z1*Z1-Z2*Z2);
      IF (ABS(AK*XMU-1.0)<1.0E-06) then XMU:=XMU*(1.0-1.0E-06);
      G:=(Z1+Z2-AK)/(Z1+Z2+AK);
      R6:=(1.0-PIZ*(B1/3.0)*(B1/3.0))*TAUT;
      E:=EXP(-AK*R6);
      H:=EXP(-R6/XMU);
      R7:=AK*AK-1.0/XMU/XMU;
      A:=1.0*(Z1*Z3+Z2*Z4-Z3/XMU)/R7;
      B:=1.0*(Z1*Z4+Z2*Z3+Z4/XMU)/R7;
      S7:=E*G;
      R8:=1.0/E;
      R9:=H*(-A);
      Q7:=S7*G-R8;
      X1:=(R9*G-(-B)*R8)/Q7;
      X2:=((-B)*S7-R9)/Q7;
      AL4:=E*X1+G*X2/E+B*H;
      TRANS:=(AL4+XMU*H)/(XMU);
END;

procedure SCATRA (WL, TAER, TRAY, GA, PIZA, XMUS, XMUV:Double;
               var DDIRTT, DDIFTT, UDIRTT, UDIFTT, SPHALBT,
                   DDIRTR, DDIFTR, UDIRTR, UDIFTR, SPHALBR,
                   DDIRTA, DDIFTA, UDIRTA, UDIFTA, SPHALBA:Double);
var
  IT : Integer;
  TT,TR,BETA1,DEL,EX1,BB,BT : Double;
  B1,PIZ,TAUT : Double;

label  997,998;

(*C     COMPUTATIONS OF THE DIRECT AND DIFFUSE TRANSMITTANCES             *)
(*C     FOR DOWNWARD AND UPWARD PATHS   AND SPHERICAL ALBEDO              *)
BEGIN
      DDIRTT:=0.;
      DDIFTT:=0.;
      UDIRTT:=0.;
      UDIFTT:=0.;
      SPHALBT:=0.;
      DDIRTR:=0.;
      DDIFTR:=0.;
      UDIRTR:=0.;
      UDIFTR:=0.;
      SPHALBR:=0.;
      DDIRTA:=0.;
      DDIFTA:=0.;
      UDIRTA:=0.;
      UDIFTA:=0.;
      SPHALBA:=0.;
 
   for IT:=1 to 3 do begin
      TT:=TAER;
      TR:=TRAY;
      IF(IT=1) then TT:=0.;
      IF(IT=2) then TR:=0.;
      IF(IT=2) AND (TAER=0.0) then GOTO 998;
      BETA1:=3.*GA;
(* Globals
                     B1, PIZ, TAUT       : Double;    *)
      TAUT:=TR+TT;
      DEL:=TR/TAUT;
      PIZ:=DEL+(1.-DEL)*PIZA;
      B1:=TT*BETA1/TAUT;
 
(*C     DOWNWARD TRANSMITTANCE                                            *)
 
      EX1:=-(TAUT)/XMUS;
      DDIRTT:=0.;
      IF(ABS(EX1)<10.0) then DDIRTT:=EXP_D(EX1);
      DDIFTT:=TRANS(XMUS,B1,PIZ,TAUT)-DDIRTT;
 
(*C     UPWARD TRANSMITTANCE                                              *)
 
      EX1:=-TAUT/XMUV;
      UDIRTT:=0.;
      IF(ABS(EX1)<10.0) then UDIRTT:=EXP(EX1);
      UDIFTT:=TRANS(XMUV,B1,PIZ,TAUT)-UDIRTT;
 
(*C     SPHERICAL ALBEDO                                                  *)
 
      BB:=0.5*(1.0-GA);
      BT:=TR+2.0*PIZA*BB*TT;
      SPHALBT:=BT*(1.+EXP_D(-2.0*BT))/(2.0+BT);
      IF (IT=2) then GOTO 997;
      IF (IT=3) then GOTO 998;
      DDIRTR:=DDIRTT;
      DDIFTR:=DDIFTT;
      UDIRTR:=UDIRTT;
      UDIFTR:=UDIFTT;
      SPHALBR:=SPHALBT;
      GOTO 998;
997 : DDIRTA:=DDIRTT;
      DDIFTA:=DDIFTT;
      UDIRTA:=UDIRTT;
      UDIFTA:=UDIFTT;
      SPHALBA:=SPHALBT;
998:
   end;

END;

var
  GA,PIZA,PFA,RORAYL,ROAERO,TRAY,TAER     : Double;
  DDIRTT, DDIFTT, UDIRTT, UDIFTT, SPHALBT : Double;
  DDIRTR, DDIFTR, UDIRTR, UDIFTR, SPHALBR : Double;
  DDIRTA, DDIFTA, UDIRTA, UDIFTA, SPHALBA : Double;
  WL : Double;
  IL : integer;

BEGIN      (* DISCOM (IAER,XMUD,XMUS,XMUV,TAER55); *)

    for IL:=1 to 10 do begin

      WL:=WLDIS^[IL];
 
(*C     COMPUTATION OF RAYLEIGH OPTICAL DEPTH AT WL                       *)
 
      ODRAYL(WL,TRAY);
(*      writeln(IL,' ',WL:7:4,' ',TRAY:7:4);   *)
      TRAYL^[IL]:=TRAY;
 
(*C     COMPUTATION OF AEROSOL OPTICAL PROPERTIES AT WL                   *)
 
      TAER:=TAER55*_EXT^[IL]/_EXT^[4];
      GA:=GASYM^[IL];
      PIZA:=OME^[IL];
      PFA:=PHASE^[IL];
(*
      writeln(IL,' ',TAER55:7:4,' ',TAER:7:4,' ',GA:7:4,' ',PIZA:7:4,' ',PFA:7:4);
*)
(*C     COMPUTATION OF ATMOSPHERIC REFLECTANCES                           *)
(*C               RORAYL IS RAYLEIGH REF                                  *)
(*C               ROAERO IS AEROSOL REF                                   *)
 
      ATMREF(IAER,XMUD,XMUS,XMUV,TRAY,TAER,GA,PIZA,PFA,
                  RORAYL,ROAERO);
(*      writeln(IL,' ',RORAYL:7:4,' ',ROAERO:7:4);
      str(WL:7:4,strw); _write(' WL : '+Strw);
      str(RORAYL:7:4,strw); _write(' RORAYL : '+Strw);
      str(ROAERO:7:4,strw); _writeln(' ROAERO : '+Strw);
*)
(*C     COMPUTATION OF SCATTERING TRANSMITANCES (DIRECT AND DIFFUSE)      *)
(*C     FIRST TIME FOR RAYLEIGH  NEXT TOTAL (RAYLEIGH+AEROSOLS)           *)
 
      SCATRA (WL, TAER, TRAY, GA, PIZA, XMUS, XMUV,
                   DDIRTT, DDIFTT, UDIRTT, UDIFTT, SPHALBT,
                   DDIRTR, DDIFTR, UDIRTR, UDIFTR, SPHALBR,
                   DDIRTA, DDIFTA, UDIRTA, UDIFTA, SPHALBA);
(*
      writeln(IL,' ','TT',DDIRTT:7:4,' ',DDIFTT:7:4);
      writeln(IL,' ','TR',DDIRTR:7:4,' ',DDIFTR:7:4);
      writeln(IL,' ','TA',DDIRTA:7:4,' ',DDIFTA:7:4);
*)
                                                  (* La =  Atmospheric Reflectance *)
      ROATM^[1,IL]:=RORAYL;        (* Rayleigh *)
      ROATM^[2,IL]:=ROAERO+RORAYL; (* Total *)
      ROATM^[3,IL]:=ROAERO;        (* Aerosol *)
                                                   (* Transmittance *)
                                        (* Downward *)
      DTDIR^[1,IL]:=DDIRTR;   (* Direct *)               (* 1. Continental *)
      DTDIF^[1,IL]:=DDIFTR;   (* Diffuse *)
      DTDIR^[2,IL]:=DDIRTT;                              (* 2. Maritime *)
      DTDIF^[2,IL]:=DDIFTT;
      DTDIR^[3,IL]:=DDIRTA;                              (* 3. Urban *)
      DTDIF^[3,IL]:=DDIFTA;
                                        (* Upward *)
      UTDIR^[1,IL]:=UDIRTR;
      UTDIF^[1,IL]:=UDIFTR;
      UTDIR^[2,IL]:=UDIRTT;
      UTDIF^[2,IL]:=UDIFTT;
      UTDIR^[3,IL]:=UDIRTA;
      UTDIF^[3,IL]:=UDIFTA;
                                                    (* Spherical Albedo *)
      SPHAL^[1,IL]:=SPHALBR;       (* Rayleigh *)
      SPHAL^[2,IL]:=SPHALBT;       (* Total *)
      SPHAL^[3,IL]:=SPHALBA;       (* Aeorosol *)
   end;

END;

procedure READ_SOLIRR(var f_spec:text);
(*C     SI (IN W/M2/MICRON) CONTAINS THE VALUES OF THE SOLAR              *)
(*C     IRRADIANCE BETWEEN 0.25 AND 4.0 MICRONS  BY STEP OF 0.005 M.      *)
(*C     THE VALUE OF THE TOTAL IRRADIANCE IS 1372 W/M2 ON THE WHOLE       *)
(*C     SPECTRUM (1358 W/M2 BETWEEN 0.25 AND 4.0 MICRONS)                 *)
var   j,k : integer;
  erreur       : Integer;

begin
      for j := 1 to 107 do begin
       readln(f_spec,S); for k := 1 to 7 do begin
          _nextpar(strw,S); val(strw, SI^[(j-1)*7+k], erreur);
       end;
      end;

      j := 107; readln(f_spec,S); for k := 1 to 2 do begin
        _nextpar(strw,S); val(strw, SI^[j*7+k], erreur);
(*        writeln(j*7+k,' ',SI^[j*7+k]);    *)
      end;
end;

procedure SOLIRR(WL:Double;var SWL:Double);
(*C     SI (IN W/M2/MICRON) CONTAINS THE VALUES OF THE SOLAR              *)
(*C     IRRADIANCE BETWEEN 0.25 AND 4.0 MICRONS  BY STEP OF 0.005 M.      *)
(*C     THE VALUE OF THE TOTAL IRRADIANCE IS 1372 W/M2 ON THE WHOLE       *)
(*C     SPECTRUM (1358 W/M2 BETWEEN 0.25 AND 4.0 MICRONS)                 *)
var
  IWL : Integer;
  PAS : Double;
begin
      PAS:=0.005;
      IWL:=round(0.5+(WL-0.250)/PAS)+1;
      IF(IWL>0) then SWL:=SI^[IWL] else SWL:=SI^[1];
end;

procedure JINV;
(*C                                                                       *)
(*C    1ST SPECTRAL BAND OF AVHRR (NOAA 9)                                *)
(*C                                                                       *)
var   i,j,k : integer;

begin
end;

procedure ODA550(IAER:integer; V:Double; var TAER55:Double);
var
      AZ, AZ1, BZ, BZ1, BNZ, BNZ1, DZ : Double;
      EV,BN5,BN51,BN23,BN231 : Double;
      i,j,k : integer;

(*C    AEROSOL OPTICAL DEPTH AT WL:=550NM                                  *)
(*C     VERTICAL REPARTITION OF AEROSOL DENSITY FOR V:=23KM                *)
(*C                     ( IN NBR OF PART/CM3 )                            *)

const
      AN23:array[1..34] of double = (
       2.828E+03, 1.244E+03, 5.371E+02, 2.256E+02, 1.192E+02,
       8.987E+01, 6.337E+01, 5.890E+01, 6.069E+01, 5.818E+01, 5.675E+01,
       5.317E+01, 5.585E+01, 5.156E+01, 5.048E+01, 4.744E+01, 4.511E+01,
       4.458E+01, 4.314E+01, 3.634E+01, 2.667E+01, 1.933E+01, 1.455E+01,
       1.113E+01, 8.826E+00, 7.429E+00, 2.238E+00, 5.890E-01, 1.550E-01,
       4.082E-02, 1.078E-02, 5.550E-05, 1.969E-08, 0.000E+00);
 
(*C     VERTICAL REPARTITION OF AEROSOL DENSITY FOR V:=5KM                 *)
(*C                     ( IN NBR OF PART/CM3 )                            *)
 
      AN5:array[1..34] of double = (
       1.378E+04, 5.030E+03, 1.844E+03, 6.731E+02, 2.453E+02,
       8.987E+01, 6.337E+01, 5.890E+01, 6.069E+01, 5.818E+01, 5.675E+01,
       5.317E+01, 5.585E+01, 5.156E+01, 5.048E+01, 4.744E+01, 4.511E+01,
       4.458E+01, 4.314E+01, 3.634E+01, 2.667E+01, 1.933E+01, 1.455E+01,
       1.113E+01, 8.826E+00, 7.429E+00, 2.238E+00, 5.890E-01, 1.550E-01,
       4.082E-02, 1.078E-02, 5.550E-05, 1.969E-08, 0.000E+00);

label 999;
 
begin

  TAER55:=0.;

  IF(V=0.0) then goto 999;
  IF(IAER=0) then goto 999;
 
  for K:=1 to 32 do begin
      DZ:=Z^[K+1]-Z^[K];
      BN5:=AN5[K];
      BN51:=AN5[K+1];
      BN23:=AN23[K];
      BN231:=AN23[K+1];
      AZ:=(115.0/18.0)*(BN5-BN23);
      AZ1:=(115.0/18.0)*(BN51-BN231);
      BZ:=(5.0*BN5/18.0)-(23.0*BN23/18.0);
      BZ1:=(5.0*BN51/18.0)-(23.0*BN231/18.0);
      BNZ:=(AZ/V)-BZ;
      BNZ1:=(AZ1/V)-BZ1;
      EV:=DZ*EXP_D((LN(BNZ)+LN(BNZ1))*0.5);
      TAER55:=TAER55 + (EV*_SIGMA*1.0E-03);
   end;

999:END;

Procedure READ_AEROSO(var f_spec:text);
var
  i,j,k : integer;
  erreur       : Integer;

procedure DUST;
(*C                                                                       *)
(*C    MODEL: DUST-LIKE                                                   *)
(*C                                                                       *)
var   i,j,k : integer;

begin
    for k := 1 to 10 do begin
      for j := 1 to 10 do begin
       readln(f_spec,S); for i := 1 to 8 do begin
          _nextpar(strw,S); val(strw, PH_D^[k,(j-1)*8+i], erreur);
       end;
      end;
      j := 10; readln(f_spec,S); for i := 1 to 3 do begin
        _nextpar(strw,S); val(strw, PH_D^[k,j*8+i], erreur);
      end;
    end;
end;

procedure WATER;
(*C                                                                       *)
(*C    MODEL: WATER-SOLUBLE                                               *)
(*C                                                                       *)
var   i,j,k : integer;
begin
    for k := 1 to 10 do begin
      for j := 1 to 10 do begin
       readln(f_spec,S); for i := 1 to 8 do begin
          _nextpar(strw,S); val(strw, PH_W^[k,(j-1)*8+i], erreur);
       end;
      end;
      j := 10; readln(f_spec,S); for i := 1 to 3 do begin
        _nextpar(strw,S); val(strw, PH_W^[k,j*8+i], erreur);
      end;
    end;
end;

procedure OCEANIC;
(*C                                                                       *)
(*C    MODEL: OCEANIC                                                     *)
(*C                                                                       *)
var   i,j,k : integer;
begin
    for k := 1 to 10 do begin
      for j := 1 to 10 do begin
       readln(f_spec,S); for i := 1 to 8 do begin
          _nextpar(strw,S); val(strw, PH_O^[k,(j-1)*8+i], erreur);
       end;
      end;
      j := 10; readln(f_spec,S); for i := 1 to 3 do begin
        _nextpar(strw,S); val(strw, PH_O^[k,j*8+i], erreur);
      end;
    end;
end;

procedure SOOT;
(*C                                                                       *)
(*C    MODEL: SOOT                                                        *)
(*C                                                                       *)
var   i,j,k : integer;
begin
    for k := 1 to 10 do begin
      for j := 1 to 10 do begin
       readln(f_spec,S); for i := 1 to 8 do begin
          _nextpar(strw,S); val(strw, PH_S^[k,(j-1)*8+i], erreur);
       end;
      end;
      j := 10; readln(f_spec,S); for i := 1 to 3 do begin
        _nextpar(strw,S); val(strw, PH_S^[k,j*8+i], erreur);
(*       writeln(PH^[k,j*8+i]);   *)
      end;
    end;
end;

begin
(*
     I: 1:=CONTINENTAL 2:=MARITIME 3:=URBAN  4:=USER
     J: 1:=DUST-LIKE 2:=WATER-SOLUBLE 3:=OCEANIC 4:=SOOT
*)

      for k := 1 to 4 do begin
        for j := 1 to 5 do begin
          readln(f_spec,S); for i := 1 to 2 do begin
            _nextpar(strw,S); val(strw, EX^[k,(j-1)*2+i], erreur);
            _nextpar(strw,S); val(strw, SC^[k,(j-1)*2+i], erreur);
          end;
        end;
      end;

      for j := 1 to 4 do begin
       readln(f_spec,S); for i := 1 to 10 do begin
          _nextpar(strw,S); val(strw, ASY^[j,i], erreur);
       end;
      end;

      for j := 1 to 10 do begin
       readln(f_spec,S); for i := 1 to 8 do begin
          _nextpar(strw,S); val(strw, C_^[(j-1)*8+i], erreur);
       end;
      end;
      j := 10; readln(f_spec,S); for i := 1 to 3 do begin
        _nextpar(strw,S); val(strw, C_^[j*8+i], erreur);
      end;
      DUST;
      WATER;
      OCEANIC;
      SOOT;

end;

Procedure AEROSO (IAER:integer;var CO_IN; var XMUD:double);

var
  CO                        : array[1..4] of double absolute CO_IN;
  NIS,SUMNI                 : double;
  _COEF,SIGM,SCA            : double;
  VI,NI,CIJ,CI              : array[1..4] of double;
  DD                        : pEX_ARR;

  i,j,k,J1,J2,L : integer;

label 998,999;

BEGIN           (* AEROSO *)

  getmem(DD,sizeof(EX_ARR));

(*
     OPTICAL PROPERTIES OF AEROSOL MODEL COMPUTED FROM SRA BASIC COMP
*)
    for L:= 1 to 10 do begin
      _EXT^[L]:=0.;
      IF(L=4) AND (IAER=0) then _EXT^[L]:=1.;
      OME^[L]:=0.;
      GASYM^[L]:=0.;
      PHASE^[L]:=0.;
    end;

      IF(IAER=0) then goto 999;

(*
     SRA BASIC COMPONENTS FOR AEROSOL MODEL, EXTINCTION COEFFICIENTS
     ARE IN KM-1.

     DUST-LIKE := 1
     WATER-SOLUBLE := 2
     OCEANIQUE := 3
     SOOT := 4
*)
      for  J:=1 to 4 do CI[J]:=CO[J];

      VI[1] := 113.983516;    VI[2] := 113.983516E-06;
      VI[3] := 5.1444150196;  VI[4] := 59.77353425E-06;

      NI[1] := 54.734;        NI[2] := 1.86855E+06;
      NI[3] := 276.05;        NI[4] := 1.80582E+06;

      for  K:=1 to 82 do begin
       IF (XMUD>=C_^[K]) AND (XMUD<C_^[K+1]) then goto 998;
      end;
      goto 999;

 998: J1:=K;
      J2:=J1+1;
      _COEF:=(XMUD-C_^[J1])/(C_^[J2]-C_^[J1]);

      for L:=1 to 10 do DD^[1,L]:=PH_D^[L,J1]+_COEF*(PH_D^[L,J1]-PH_D^[L,J2]);
      for L:=1 to 10 do DD^[2,L]:=PH_W^[L,J1]+_COEF*(PH_W^[L,J1]-PH_W^[L,J2]);
      for L:=1 to 10 do DD^[3,L]:=PH_O^[L,J1]+_COEF*(PH_O^[L,J1]-PH_O^[L,J2]);
      for L:=1 to 10 do DD^[4,L]:=PH_S^[L,J1]+_COEF*(PH_S^[L,J1]-PH_S^[L,J2]);

(*
     SUMM OF THE CI/VI CALCULATION
*)
      SUMNI:=0.;
      SIGM:=0.;
      for I:=1 to 4 do SIGM:=SIGM+CI[I]/VI[I];

(*
     CIJ COEFFICIENTS CALCULATION
*)

      for J:=1 to 4 do begin
        CIJ[J]:=(CI[J]/VI[J]/SIGM);
        SUMNI:=SUMNI+CIJ[J]/NI[J];
      end;
      NIS:=1.E+00/SUMNI;

(*
     MIXING PARAMETERS CALCULATION
*)

    for L:=1 to 10 do begin
      SCA:=0.;
      for J:=1 to 4  do begin
        _EXT^[L]:=EX^[J,L]*CIJ[J]+_EXT^[L];
        SCA:=SC^[J,L]*CIJ[J]+SCA;
        GASYM^[L]:=SC^[J,L]*CIJ[J]*ASY^[J,L]+GASYM^[L];
        PHASE^[L]:=SC^[J,L]*CIJ[J]*DD^[J,L]+PHASE^[L];
      end;
      OME^[L]:=SCA/_EXT^[L];
      _EXT^[L]:=_EXT^[L]*NIS;
      GASYM^[L]:=GASYM^[L]/SCA;
      PHASE^[L]:=PHASE^[L]/SCA;
    end;

999:
  freemem(DD,sizeof(EX_ARR));
END;

Procedure  ATMOSPHERE_MODEL(iatmo:integer) ;
var   i,j,k : integer;
  erreur       : Integer;

  f_atm    : file;
  buf_data : array[1..5] of Double;

label 999;

BEGIN
(*
    MODEL: Mid-latitude summer MC CLATCHEY
    MODEL: US STANDARD 62 MC CLATCHEY
*)
(*
     for j := 1 to 3 do begin
       readln(f_spec,S); for i := 1 to 9 do begin
          _nextpar(strw,S); val(strw, Z^[(j-1)*9+i], erreur);
       end;
     end;
     j := 3; readln(f_spec,S); for i := 1 to 7 do begin
        _nextpar(strw,S); val(strw, Z^[j*9+i], erreur);
     end;

     for j := 1 to 5 do begin
       readln(f_spec,S); for i := 1 to 6 do begin
          _nextpar(strw,S); val(strw, P^[(j-1)*6+i], erreur);
       end;
     end;
     j := 5; readln(f_spec,S); for i := 1 to 4 do begin
        _nextpar(strw,S); val(strw, P^[j*6+i], erreur);
     end;

     for j := 1 to 5 do begin
       readln(f_spec,S); for i := 1 to 6 do begin
          _nextpar(strw,S); val(strw, T^[(j-1)*6+i], erreur);
       end;
     end;
     j := 5; readln(f_spec,S); for i := 1 to 4 do begin
        _nextpar(strw,S); val(strw, T^[j*6+i], erreur);
     end;

     for j := 1 to 5 do begin
       readln(f_spec,S); for i := 1 to 6 do begin
          _nextpar(strw,S); val(strw, WH^[(j-1)*6+i], erreur);
       end;
     end;
     j := 5; readln(f_spec,S); for i := 1 to 4 do begin
        _nextpar(strw,S); val(strw, WH^[j*6+i], erreur);
     end;

     for j := 1 to 5 do begin
       readln(f_spec,S); for i := 1 to 6 do begin
          _nextpar(strw,S); val(strw, WO^[(j-1)*6+i], erreur);
       end;
     end;
     j := 5; readln(f_spec,S); for i := 1 to 4 do begin
        _nextpar(strw,S); val(strw, WO^[j*6+i], erreur);
     end;

{JMM} assign(f_atm,ExePath+'ATMO6.DAT');
{JMM} rewrite(f_atm,1);
      ier := ioresult;
      if ier <> 0 then goto 999;

      for j := 1 to 34 do begin
        buf_data[1] := Z^[j];  buf_data[2] := P^[j];  buf_data[3] := T^[j];
        buf_data[4] := WH^[j]; buf_data[5] := WO^[j];
        blockwrite(f_atm,buf_data,40); { Double = 8 bytes }
      end;

{JMM} close(f_atm);
{$I+}
*)
{$I-}
      str(iatmo:1,strw);
{JMM} assign(f_atm,ExePath+'ATMO'+strw+'.DAT');
{JMM} reset(f_atm,1);
      ier := ioresult;
      if ier <> 0 then goto 999;

      for j := 1 to 34 do begin
        blockread(f_atm,buf_data,40,ier); { Double = 8 bytes }
        Z^[j] := buf_data[1]; P^[j] := buf_data[2]; T^[j] := buf_data[3];
        WH^[j] := buf_data[4]; WO^[j] := buf_data[5];
      end;

{JMM} close(f_atm);
{$I+}

999:END;

Function ASIN(X:Double):Double;
begin ASIN := ArcTan (x/sqrt (1-sqr (x))); end;

Function ACOS(X:Double):Double;
begin ACOS := ArcTan (sqrt (1-sqr (x)) /x); end;


Function POSSOL (MONTH,JDAY:integer;TU,XLON,XLAT:double;
                        var ASOL,PHI0:double):integer;
(*
     SOLAR POSITION (ZENITHAL ANGLE ASOL,AZIMUTHAL ANGLE PHI0
                     IN DEGREES)
     JDAY IS THE NUMBER OF THE DAY IN THE MONTH
*)
label 999;
var
  PI,FAC,PI2            : Double;
  J                     : integer;
  TSM,XLO,XLA,XJ        : Double;
  A1,A2,ET              : Double;
  TSV,AH                : Double;
  A3,DELTA_,AMUZERO     : Double;
  ELEV,AZ,CAZ,AZIM      : Double;

  begin
    POSSOL := -1;

      PI  :=3.1415926;
      FAC :=PI/180.;

      IF (MONTH <= 2)then J:=31*(MONTH-1)+JDAY;
      IF  MONTH in [3..8] then J:=31*(MONTH-1)-((MONTH-1)div 2)-2+JDAY;
      IF (MONTH > 8) then J:=31*(MONTH-1)-((MONTH-2) div 2)-2+JDAY;
(*
    MEAN SOLAR TIME
*)
      TSM:=TU+XLON/15.;
      XLO:=XLON*FAC;
      XLA:=XLAT*FAC;
      XJ:=J;
(*
    TIME EQUATION (IN mn.dec)
*)
      A1:=(1.00554*XJ-6.28306)*FAC;
      A2:=(1.93946*XJ+23.35089)*FAC;
      ET:=-7.67825*SIN(A1)-10.09176*SIN(A2);
(*
      TRUE SOLAR TIME
*)
      TSV:=TSM+ET/60.;
      TSV:=TSV-12.;
(*
      HOUR ANGLE
*)
      AH:=TSV*15.*FAC;
(*
      SOLAR DECLINATION   (IN RADIAN)
*)
      A3:=(0.9683*XJ-78.00878)*FAC;
      DELTA_:=23.4856*SIN(A3)*FAC;
(*
      ELEVATION,AZIMUTH
*)
      AMUZERO:=SIN(XLA)*SIN(DELTA_)+COS(XLA)*COS(DELTA_)*COS(AH);
      ELEV:=ASIN(AMUZERO);
      AZ:=COS(DELTA_)*SIN(AH)/COS(ELEV);
      CAZ:=(-COS(XLA)*SIN(DELTA_)+SIN(XLA)*COS(DELTA_)*COS(AH))/COS(ELEV);
      AZIM:=ASIN(AZ);
      IF CAZ<=0.0 then AZIM:=PI-AZIM;
      IF (CAZ>0.0) and (AZ<=0.0) then AZIM:=2*PI+AZIM ;
      AZIM:=AZIM+PI;
      PI2:=2*PI;
      IF AZIM>PI2 then AZIM:=AZIM-PI2;
      ELEV:=ELEV*180./PI;
(*
      CONVERSION IN DEGREES
*)
      ASOL:=90.-ELEV;
      PHI0:=AZIM/FAC;
      IF(ASOL>=90) then begin
        Loadstring(HInstance, stW_SpecRIS+Lang, MessageText, 79);
        Loadstring(HInstance, stW_SpecTIT+Lang, TitleText, 79);
        MessageBox(0, MessageText, TitleText, mb_IconStop);
        goto 999;
      end;

    POSSOL := 0;
999:end;

Function  POSMTO (MONTH,JDAY:integer;TU:Double;NC,NL:integer;
                        var ASOL,PHI0,AVIS,PHIV,XLON,XLAT:double):integer;
var
  YR,XR,ALTI,RE,AAA,RP,PI,CDR,CRD,DELTAX,DELTAY : Double;
  X,Y,RS,TANX,TANY,VAL1,VAL2,YK,COSX2  : Double;
  SN,ZT,XT,YT,TETA,YLAT,YLON,GAM : Double;

label 999;

begin

   POSMTO := -1;

(*C     METEOSAT DEFINITION                                               *)
 
      YR:=NL-1250.5;
      XR:=NC-2500.5;
      ALTI:=42164.0-6378.155;
      RE:=6378.155;
      AAA:=1./297.;
      RP:=RE/(1.+AAA);
      PI:=3.1415926;
      CDR:=PI/180.;
      CRD:=180./PI;
      DELTAX:=18.0/5000.0;
      DELTAY:=18.0/2500.0;
      X:=XR*DELTAX*CDR;
      Y:=YR*DELTAY*CDR;
      RS:=RE+ALTI;
      TANX:=SIN(X)/COS(X);
      TANY:=SIN(Y)/COS(Y);
      VAL1:=1.0+(TANX*TANX);
      VAL2:=1.0+(TANY*(1.0+AAA))*(TANY*(1.0+AAA));
      YK:=RS/RE;
      COSX2:=1./(VAL1*VAL2);
      IF (1./COSX2) > ((YK*YK)/(YK*YK-1.0)) then begin
        Loadstring(HInstance, stW_SpecLAT+Lang, MessageText, 79);
        Loadstring(HInstance, stW_SpecTIT+Lang, TitleText, 79);
        MessageBox(0, MessageText, TitleText, mb_IconStop);
        goto 999;
      end;
      SN:=(RS-(RE*(SQRT((YK*YK)-(YK*YK-1.0)*(1.0/COSX2)))))/(1.0/COSX2);
      ZT:=RS-SN;
      XT:=-(SN*TANX);
      YT:=SN*TANY/COS(X);
      TETA:=ASIN(YT/RP);
      YLAT:=(ARCTAN(((SIN(TETA)/COS(TETA))*RP)/RE));
      YLON:=ARCTAN(XT/ZT);

      XLAT:=YLAT*CRD;
      XLON:=YLON*CRD;

      ier := POSSOL(MONTH,JDAY,TU,XLON,XLAT,ASOL,PHI0);
      if ier <> 0 then goto 999;
 
      YLON:=XLON*PI/180.;
      YLAT:=XLAT*PI/180.;
      GAM:=SQRT(((1.0/COSX2)-1.0)*COSX2);
      AVIS:=ASIN((1.0+ALTI/RE)*(GAM));
      AVIS:=AVIS*180./PI;
      PHIV:=ARCTAN(SIN(YLON)/(COS(YLON)*SIN(YLAT)))+PI;
      PHIV:=PHIV*180./PI;

   POSMTO := 0;
999: END;

Function  POSGE (MONTH,JDAY:integer;TU:Double;NC,NL:integer;
                        var ASOL,PHI0,AVIS,PHIV,XLON,XLAT:double):integer;
var
  YR,XR,ALTI,RE,AAA,RP,PI,CDR,CRD,DELTAX,DELTAY : Double;
  X,Y,RS,TANX,TANY,VAL1,VAL2,YK,COSX2  : Double;
  SN,ZT,XT,YT,TETA,YLAT,YLON,GAM : Double;

label 999;

begin
  POSGE := -1;
(*C     GOES EAST DEFINITION                                              *)
 
      YR:=NL-8665.5;
      XR:=NC-6498.5;
      ALTI:=42107.0-6378.155;
      RE:=6378.155;
      AAA:=1./297.;
      RP:=RE/(1.+AAA);
      PI:=3.1415926;
      CDR:=PI/180.;
      CRD:=180./PI;
      DELTAX:=18.0/12997.0;
      DELTAY:=20.0/17331.0;
      X:=XR*DELTAX*CDR;
      Y:=YR*DELTAY*CDR;
      RS:=RE+ALTI;
      TANX:=SIN(X)/COS(X);
      TANY:=SIN(Y)/COS(Y);
      VAL1:=1.0+(TANX*TANX);
      VAL2:=1.0+(TANY*(1.0+AAA))*(TANY*(1.0+AAA));
      YK:=RS/RE;
      COSX2:=1./(VAL1*VAL2);
      IF (1./COSX2) > ((YK*YK)/(YK*YK-1.0)) then begin
        Loadstring(HInstance, stW_SpecLAT+Lang, MessageText, 79);
        Loadstring(HInstance, stW_SpecTIT+Lang, TitleText, 79);
        MessageBox(0, MessageText, TitleText, mb_IconStop);
        goto 999;
      end;
      SN:=(RS-(RE*(SQRT((YK*YK)-(YK*YK-1.0)*(1.0/COSX2)))))/(1.0/COSX2);
      ZT:=RS-SN;
      XT:=-(SN*TANX);
      YT:=SN*TANY/COS(X);
      TETA:=ASIN(YT/RP);
      YLAT:=(ARCTAN(((SIN(TETA)/COS(TETA))*RP)/RE));
      YLON:=ARCTAN(XT/ZT);

      XLAT:=YLAT*CRD;
      XLON:=YLON*CRD-75.;
 
      ier := POSSOL(MONTH,JDAY,TU,XLON,XLAT,ASOL,PHI0);
      if ier <> 0 then goto 999;

      YLON:=XLON*PI/180.+75.*CDR;
      YLAT:=XLAT*PI/180.;
      GAM:=SQRT(((1./COSX2)-1.0)*COSX2);
      AVIS:=ASIN((1.+ALTI/RE)*(GAM));
      AVIS:=AVIS*180./PI;
      PHIV:=ARCTAN(SIN(YLON)/(COS(YLON)*SIN(YLAT)))+PI;
      PHIV:=PHIV*180./PI;

    POSGE := 0;
999:END;

Function  POSGW (MONTH,JDAY:integer;TU:Double;NC,NL:integer;
                        var ASOL,PHI0,AVIS,PHIV,XLON,XLAT:double):integer;
var
  YR,XR,ALTI,RE,AAA,RP,PI,CDR,CRD,DELTAX,DELTAY : Double;
  X,Y,RS,TANX,TANY,VAL1,VAL2,YK,COSX2  : Double;
  SN,ZT,XT,YT,TETA,YLAT,YLON,GAM : Double;

label 999;

begin
    POSGW := -1;

(*C     GOES WEST DEFINITION                                              *)
 
      YR:=NL-8665.5;
      XR:=NC-6498.5;
      ALTI:=42147.0-6378.155;
      RE:=6378.155;
      AAA:=1.0/297.;
      RP:=RE/(1.+AAA);
      PI:=3.1415926;
      CDR:=PI/180.;
      CRD:=180./PI;
      DELTAX:=18.0/12997.0;
      DELTAY:=20.0/17331.0;
      X:=XR*DELTAX*CDR;
      Y:=YR*DELTAY*CDR;
      RS:=RE+ALTI;
      TANX:=SIN(X)/COS(X);
      TANY:=SIN(Y)/COS(Y);
      VAL1:=1.0+(TANX*TANX);
      VAL2:=1.0+(TANY*(1.0+AAA))*(TANY*(1.0+AAA));
      YK:=RS/RE;
      COSX2:=1./(VAL1*VAL2);
      IF (1./COSX2) > ((YK*YK)/(YK*YK-1.0)) then begin
        Loadstring(HInstance, stW_SpecLAT+Lang, MessageText, 79);
        Loadstring(HInstance, stW_SpecTIT+Lang, TitleText, 79);
        MessageBox(0, MessageText, TitleText, mb_IconStop);
        goto 999;
      end;
      SN:=(RS-(RE*(SQRT((YK*YK)-(YK*YK-1.0)*(1.0/COSX2)))))/(1.0/COSX2);
      ZT:=RS-SN;
      XT:=-(SN*TANX);
      YT:=SN*SIN(Y)/COS(Y)*COS(X);
      TETA:=ASIN(YT/RP);
      YLAT:=(ARCTAN(((SIN(TETA)/COS(TETA))*RP)/RE));
      YLON:=ARCTAN(XT/ZT);

      XLAT:=YLAT*CRD;
      XLON:=YLON*CRD-135.;
 
      ier := POSSOL(MONTH,JDAY,TU,XLON,XLAT,ASOL,PHI0);
      if ier <> 0 then goto 999;

      YLON:=XLON*PI/180.+135.*CDR;
      YLAT:=XLAT*PI/180.;
      GAM:=SQRT(((1./COSX2)-1.0)*COSX2);
      AVIS:=ASIN((1.+ALTI/RE)*(GAM));
      AVIS:=AVIS*180./PI;
      PHIV:=ARCTAN(SIN(YLON)/(COS(YLON)*SIN(YLAT)))+PI;
      PHIV:=PHIV*180./PI;

    POSGW := 0;
999:END;

Function  POSNOA (MONTH,JDAY:integer;TU:Double;NC:integer;XLONAN,HNA:Double;
                        var ASOL,PHI0,AVIS,PHIV,XLON,XLAT:double):integer;
var
  PI,R,AI,AN,YLONAN,T,HNAM,U,DELT,D,X,Y,Z : Double;
  YLAT,COSY,SINY,YLON,YLO1 : Double;
  ZLAT,ZLON,XNUM,XDEN : Double;

label 10,999;

begin
  POSNOA := -1;

(*C     NOAA 6 DEFINITION                                                 *)
(*C     ORBITE INCLINATION AI IN RADIANS                                  *)
(*C     HOR MOUV IN RAD/S  AN                                             *)
(*C     H/R:=860/6378                                                      *)
 
      PI:=3.1415926;
      R:=860./6378.155;
      AI:=98.96*PI/180.;
      AN:=360.*PI/(6119.*180.0);
      YLONAN:=XLONAN*PI/180.0;
      T:=TU*3600;
      HNAM:=HNA;
      HNAM:=HNAM*3600;
      U:=T-HNAM;
      U:=U*AN;
      DELT:=((NC-1024)*55.385/1024.0);
      DELT:=DELT*PI/180.;
      AVIS:=ASIN((1+R)*SIN(DELT));
      D:=AVIS-DELT;
      X:=COS(D)*COS(U);
      Y:=COS(D)*COS(AI)*SIN(U)-SIN(AI)*SIN(D);
      Z:=COS(D)*SIN(AI)*SIN(U)+COS(AI)*SIN(D);
      YLAT:=ASIN(Z);
      COSY:=COS(D)*COS(U)/COS(YLAT);
      SINY:=Y/COS(YLAT);
      YLON:=ASIN(SINY);
      IF(COSY > 0.0) then GOTO 10;
      IF(SINY>0.0) then YLON:=PI-YLON;
      IF(SINY<=0.0) then YLON:=-(PI+YLON);
 10:  YLO1:=YLON+YLONAN-(T-HNAM)*2.*PI/86400.;
      XLAT:=YLAT*180./PI;
      XLON:=YLO1*180./PI;
 
      ier := POSSOL(MONTH,JDAY,TU,XLON,XLAT,ASOL,PHI0);
      if ier <> 0 then goto 999;

      ZLAT:=ASIN(SIN(AI)*SIN(U));
(*      ZLON:=ATAN2(COS(AI)*SIN(U) COS(U));    *)
      XNUM:=SIN(ZLON-YLON)*COS(ZLAT)/SIN(ABS(D));
      XDEN:=(SIN(ZLAT)-SIN(YLAT)*COS(D))/COS(YLAT)/SIN(ABS(D));
      PHIV:=0.;
(*      IF(NC<>1024)PHIV:=ATAN2(XNUM XDEN);  *)
      PHIV:=PHIV*180./PI;
      AVIS:=ABS(AVIS);
      AVIS:=AVIS*180./PI;

   POSNOA := 0;
999: END;


Function POSPO (MONTH,JDAY:integer;TU,XLON,XLAT:double;
                        var ASOL,PHI0,AVIS,PHIV:double):integer;
(*C     SPOT DEFINITION                                                   *)
(*C     WARNING !!!                                                       *)
(*C     XLON AND XLAT ARE THE COORDINATES OF THE SCENE CENTER.            *)
 
begin
      AVIS:=0.;
      PHIV:=0.;

      ier := POSSOL(MONTH,JDAY,TU,XLON,XLAT,ASOL,PHI0);
      POSPO := ier;
end;

Function POSLAN (MONTH,JDAY:integer;TU,XLON,XLAT:double;
                        var ASOL,PHI0,AVIS,PHIV:double):integer;

(*
     LANDSAT5 DEFINITION
     WARNING !!!
     XLON AND XLAT ARE THE COORDINATES OF THE SCENE CENTER.
*)

begin
      AVIS:=0.;
      PHIV:=0.;

      ier := POSSOL(MONTH,JDAY,TU,XLON,XLAT,ASOL,PHI0);
      POSLAN := ier;
end;

procedure READ_WLDISC(var f_spec:text);
var   j,k : integer;
  erreur       : Integer;

begin
      for j := 1 to 2 do begin
       readln(f_spec,S); for k := 1 to 5 do begin
          _nextpar(strw,S); val(strw, WLDISC^[(j-1)*5+k], erreur);
       end;
      end;
end;

procedure READ_CO3(var f_spec:text);
var   j,k : integer;
  erreur       : Integer;

begin
      for j := 1 to 17 do begin
       readln(f_spec,S); for k := 1 to 6 do begin
          _nextpar(strw,S); val(strw, CO3^[(j-1)*6+k], erreur);
       end;
      end;
end;

procedure READ_CCH2O(var f_spec:text);
var   j,k : integer;
  erreur       : Integer;

begin
      for j := 1 to 3 do begin
       readln(f_spec,S); for k := 1 to 5 do begin
          _nextpar(strw,S); val(strw, CCH2O^[(j-1)*5+k], erreur);
       end;
      end;
end;

end.